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1  Abstract 

The  investigators  report  on  their  findings,  recent  publications  and  presentations  in  the  areas 
of  lenslet  array  imaging,  wavefront  encoding,  and  non-negative  matrix  factorization  for  ma¬ 
terial  component  (end-member)  identification.  Lenslet  arrays  enable  a  number  of  imaging 
modalities,  including  amplitude  diversity,  polarization  diversity,  wavelength  diversity  (multi- 
spectral  data)  and  phase  diversity.  Each  of  these  techniques  extends  traditional  imaging  by 
modifying  the  data  acquisition  to  implicitly  capture  features  that  would  otherwise  be  un¬ 
detectable.  Post-processing  converts  implicitly  captured  or  encoded  information  to  a  form 
suitable  for  human  or  automated  identification  tasks. 

The  problem  of  material  identification  from  multi-spectral  image  data  (blind  source  sepa¬ 
ration)  can  be  formulated  as  a  non-negative  matrix  factorization  problem  or  as  a  tensor 
factorization  problem.  Non-uniqueness  and  numerical  stability  are  frequently  difficult  issues. 
We  report  here  on  recent  progress  on  using  stochastic  constraints  for  improving  the  numerical 
stability  of  the  computation. 

Current  (2009)  and  prior  results  relating  to  pupil  phase  encoding  were  adapted  to  mitigate 
a  frequency-agile  pulsed  laser  attack  against  CCD-based  cameras.  The  investigators  per¬ 
formed  a  phase  one  SBIR  with  Agiltron  Corp.  to  explore  possible  commercialization  of  pupil 
phase  encoding  technology.  While  this  approach  was  not  selected  for  phase  two  funding,  the 
investigators  showed  that  this  approach  is  viable. 


2  Publications 

2.1  Publications  in  Peer-review  Journals 

1.  V.P.  Pauca,  K.  Smith,  A.  Ross,  and  J.  van  der  Gracht.  “Wavefront  Coded  Iris  Biometric 
Systems,”  Encyclopedia  of  Biometrics,  Li,  Stan  Z.  (Ed.),  Springer  (ISBN:  978-0-387- 
73002-8),  2009. 

2.  P.  Zhang,  H.  Wang,  R.J.  Plemmons,  and  V.P.  Pauca.  “Tensor  Methods  for  Hyperspec- 
tral  Analysis  of  Space  Objects,  A  Material  Identification  Study,”  Journal  of  the  Optical 
Society  of  America,  A.,  Vol.  25,  No.  12,  pp.  3001-3012,  2008. 
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3.  M.W.  Berry,  M.  Browne,  A.N.  Langville,  V.P.  Pauca,  and  R.  J.  Plemmons.  “Algorithms 
and  Applications  for  Approximate  Nonnegative  Matrix  Factorization,”  Computational 
Statistics  &  Data  Analysis,  Vol.  52,  pp.  155-173,  2007. 

4.  V.P.  Pauca,  J.  Piper,  and  R.J.  Plemmons.  “Nonnegative  Matrix  Factorization  for 
Spectral  Data  Analysis,”  Linear  Algebra  and  its  Applications,  vol.  416,  pp.  29-47, 
2006. 

5.  D.  Chen  and  R.  Plemmons,  “Nonnegativity  Constraints  in  Numerical  Analysis,”  To 
appear  in  the  Conference  Proceedings  of  the  Symposium  on  the  Birth  of  Numerical 
Analysis,  held  Leuven  Belgium,  October  2007.  World  Scientific  Press,  A.  Bultheel  and 
R.  Cools,  Eds. 

2.2  Non-Peer-Reviewed  Conference  Proceedings 

6.  J.  van  der  Gracht,  L.  Zhang,  T.  Torgersen,  and  P.  Pauca.  “Pupil  Phase  Encoding  for 
Mitigation  of  Laser-Induced  Saturation  in  Imaging  Sensors”.  Computational  Optical 
Sensing  and  Imaging,  OSA  Technical  Digest,  San  Jose,  California,  October  2009. 

7.  V.P.  Pauca  ,  D.  Chen,  J.  van  der  Gracht,  R.J.  Plemmons,  S.  Prasad,  and  T.C.  Torg¬ 
ersen.  “Pupil  Phase  Encoding  for  Multi-Aperture  Imaging,”  to  appear  in  Proc.  SPIE, 
Advanced  Signal  Processing  Algorithms,  Architectures,  and  Implementations  XVIII, 
San  Diego,  CA,  Aug.  2008. 

8.  K.N.  Smith,  V.P.  Pauca,  A.  Ross,  T.  Torgersen,  and  M.  C.  King.  “Extended  Evaluation 
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B.  Gray,  V.  Pauca,  T.  Torgersen,  J.  van  der  Gracht,  and  G.  Behrmann.  “PERIODIC: 
Integrated  Computational  Array  Imaging  Technology,”  OSA  Technical  Digest,  Confer¬ 
ence  on  Computational  Optical  Sensing  and  Imaging  (COSI),  Vancouver,  2007. 

10.  (extended  abstract)  R.  Barnard,  B.  Gray,  V.  Pauca,  T.  Torgersen,  M.  Mirotznik,  J.  van 
der  Gracht,  R.  Plemmons,  G.  Behrmann,  S.  Mathews,  and  S.  Prasad.  “PERIODIC: 
State-of-the-Art  Array  Imaging  Technology,”  Proc.  2007  ACM  Southeast  Conference, 
pp.  544-546,  Mar.,  2007. 

11.  (extended  abstract)  Q.  Zhang,  H.  Wang,  R.  Plemmons,  and  V.  Pauca.  “Spectral 
Unmixing  Using  Nonnegative  Tensor  Factorization,”  Proc.  2007  ACM  Southeast  Con¬ 
ference,  pp.  531-532,  Mar.,  2007. 

12.  R.  Barnard,  V.P.  Pauca,  T.C.  Torgersen,  R.J.  Plemmons,  S.  Prasad,  J.  van  der  Gracht, 
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Image  Reconstruction  from  Low-Resolution  Imagery,”  Proc.  SPIE,  Advanced  Signal 
Processing  Algorithms,  Architectures,  and  Implementations  XVI,  vol.  6313,  pp.  1-13, 
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2.3  Book  Chapter 


13.  M.  Ng  and  R.  Plemmons,  “Blind  Deconvolution  and  Structured  Matrix  Computations 
with  Applications  to  Array  Imaging”,  46  Page  Chapter  in  Blind  Image  Deconvolution: 
Theory  and  Applications,  Edited  by  P.  Campisi  and  K.  Egiazarian,  CRC  Press,  Octo¬ 
ber  2007. 


3  Presentations 

•  J.  van  der  Gracht,  Lei  Zhang,  V.  P.  Pauca,  and  T.  Torgersen,  “Pupil  Phase  Encoding 
for  Mitigation  of  Laser-induced  Saturation  in  Imaging  Sensors” ,  Computational  Optical 
Sensing  and  Imaging  (COSI),  San  Jose,  CA,  October  13-15,  2009. 

•  Technical  Conference:  October  13-15,  2009  R.  Plemmons,  (Invited),  “Nonnegativity 
Constraints  in  Numerical  Analysis”,  Symposium  on  :  The  Birth  of  Numerical  Analysis, 
Leuven,  Belgium,  October  2007. 

•  R.  Plemmons,  (Invited),  “Integrated  Computational  Lenslet  Array  Systems  Applied 
to  Bioimaging”,  Special  Semester  on  Quantitative  Biology  Analyzed  by  Mathematical 
Methods,  Linz,  Austria,  November  2007. 

•  R.  Plemmons,  (Invited),  “Matrix  and  Tensor  Methods  for  Space  Situational  Aware¬ 
ness”,  Colloquium  at  Louvain-la-Neuve  University,  Belgium,  March  2008. 

•  R.  Plemmons,  (Invited),  “Mathematical  and  Algorithms  Challenges  for  Processing  Hy- 
perspectral  Data” ,  Croucher  Advanced  Study  Institute  on  Mathematical  and  Algorith¬ 
mic  Challenges  for  Modeling  and  Analyzing  Modern  Data  Sets,  Hong  Kong,  China, 
April  2008. 

•  R.  Plemmons,  (Invited),  “Nonnegativity  Constraints  for  Linear  Algebra  in  Scientific 
Computation”,  Householder  Symposium,  Berlin,  Germany,  May  2008. 

•  V.P.  Pauca  (Invited),  “PERIODIC:  Multi- Aperture  Multi-Diversity  Imaging,”  Commu¬ 
nity  Academic  Summit,  Sponsored  by  the  Defense  Intelligence  Agency,  The  National 
Geospatial  Intelligence  Agency,  and  the  Disrupted  Technologies  Office,  Williamsburg, 
VA,  June,  2007. 

•  V.P.  Pauca.  (Invited),  “Computational  Imaging  Systems:  The  Growing  Intersection 
between  Scientific  Disciplines,”  Universidad  Andina  del  Cusco,  Cusco,  Peru,  July  2008. 

•  V.P.  Pauca.  “Pupil  Phase  Encoding  for  Multi-Aperture  Imaging,”  SPIE  Advanced 
Signal  Processing  Algorithms,  Architectures,  and  Implementations  XVIII,  San  Diego, 
CA,  August  2008. 

•  T.  Torgersen,  “PERIODIC:  Multi-Spectral  Imaging  for  Rapid  Evaluation  of  Thermal 
Injury”,  presented  at  the  2008  SIAM  Annual  Meeting,  held  July  7-11,  San  Diego,  CA. 

•  T.  Torgersen,  “PERIODIC:  Lenslet  Array  Imaging  Technology”,  presented  at  the  6th 
International  Congress  for  Industrial  and  Applied  Mathematics,  held  July  16-20,  2007, 
Zurich,  Switzerland. 
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•  T.  Torgersen,  “PERIODIC:  Digital  Superresolution  using  Lenslet  Array  Technology” 
NATO  ET051  Lab:  Industry  Day  Briefings  on  Super-Resolution  and  Advanced  Signal 
Processing,  held  September  26,  2007  at  the  US  Army  Night  Vision  Labs,  Ft.  Belvoir, 
VA. 

4  Scientific  Progress  and  Accomplishments 

4.1  Summary  of  Results 

We  report  on  our  scientific  results  in  abstract  summary  form  given  below: 

Pupil  Phase  Encoding 

Digital  super-resolution  refers  to  computational  techniques  that  exploit  the  generalized 
sampling  theorem  to  extend  image  resolution  beyond  the  pixel  spacing  of  the  detector, 
but  not  beyond  the  optical  limit  (Nyquist  spatial  frequency)  of  the  lens.  The  approach 
to  digital  super-resolution  taken  by  our  multi-lenslet  camera  project  is  to  solve  a  forward 
model  which  describes  the  effects  of  sub-pixel  shifts,  optical  blur,  and  detector  sampling 
as  a  product  of  matrix  factors.  The  associated  system  matrix  is  often  ill-conditioned, 
and  convergence  of  iterative  methods  to  solve  for  the  high-resolution  image  may  be 
slow. 

We  have  investigated  and  published1  the  use  of  pupil  phase  encoding  in  a  multi-lenslet 
camera  system  as  a  means  to  physically  precondition  and  regularize  the  computational 
super-resolution  problem.  This  is  an  integrated  optical-digital  approach  that  has  been 
previously  demonstrated  with  cubic  type  and  pseudo-random  phase  elements.  Tradi¬ 
tional  multi-frame  phase  diversity  for  imaging  through  atmospheric  turbulence  uses  a 
known  smooth  phase  perturbation  to  help  recover  a  time  series  of  point  spread  func¬ 
tions  corresponding  to  random  phase  errors.  In  the  context  of  a  multi-lenslet  camera 
system,  a  known  pseudo-random  or  cubic  phase  error  may  be  used  to  help  recover  an 
array  of  unknown  point  spread  functions  corresponding  to  manufacturing  and  focus 
variations  among  the  lenslets. 

The  use  of  pupil  phase  encoding  was  also  investigated  as  a  possible  means  to  mitigate 
a  frequency  agile  laser  attack  on  a  CCD-based  camera.  Results  are  reported  in  [6]. 
An  alternative  approach  based  on  a  lenslet  array  design  was  also  shown  effective  in 
simulation;  however  this  technical  approach  also  implied  a  significantly  reduced  field  of 
view. 

Biometric  Identification 

We  investigated  the  use  of  a  novel  multi-lens  imaging  system  in  the  context  of  biometric 
identification2,  and  more  specifically,  for  iris  recognition.  Multi-lenslet  cameras  offer  a 
number  of  significant  advantages  over  standard  single-lens  camera  systems,  including 
thin  form-factor  and  wide  angle  of  view.  By  using  appropriate  lenslet  spacing  relative 

1See  references  7  in  section  2.2. 

2See  reference  12  in  section  §2.2. 
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to  the  detector  pixel  pitch,  the  resulting  ensemble  of  images  implicitly  contains  subject 
information  at  higher  spatial  frequencies  than  those  present  in  a  single  image.  Addi¬ 
tionally,  a  multi-lenslet  approach  enables  the  use  of  observational  diversity,  including 
phase,  polarization,  neutral  density,  and  wavelength  diversities.  For  example,  post¬ 
processing  multiple  observations  taken  with  differing  neutral  density  filters  yields  an 
image  having  an  extended  dynamic  range.  Our  research  group  has  developed  several 
multi-lens  camera  prototypes  for  the  investigation  of  such  diversities. 

We  presented  techniques  for  computing  a  high-resolution  reconstructed  image  from 
an  ensemble  of  low-resolution  images  containing  sub-pixel  level  displacements.  The 
quality  of  a  reconstructed  image  is  measured  by  computing  the  Hamming  distance 
between  the  Daugman4  iris  code  of  a  conventional  reference  iris  image,  and  the  iris 
code  of  a  corresponding  reconstructed  image.  We  present  numerical  results  concerning 
the  effect  of  noise  and  defocus  blur  in  the  reconstruction  process  using  simulated  data 
and  report  preliminary  work  on  the  reconstruction  of  actual  iris  data  obtained  with 
our  camera  prototypes. 

Non-Negative  Matrix  Factorization 


In  reference  3  (section  §2.1  we  discuss  the  development  and  use  of  low-rank  approx¬ 
imate  nonnegative  matrix  factorization  (NMF)  algorithms  for  feature  extraction  and 
identification  in  the  fields  of  text  mining  and  spectral  data  analysis.  The  evolution  and 
convergence  properties  of  hybrid  methods  based  on  both  sparsity  and  smoothness  con¬ 
straints  for  the  resulting  nonnegative  matrix  factors  are  discussed.  The  interpretation 
of  NMF  outputs  in  specific  contexts  are  provided  along  with  opportunities  for  future 
work  in  the  modification  of  NMF  algorithms  for  large-scale  and  time-varying  datasets. 

Reference  5  (section  2.1  gives  a  survey  of  the  development  of  algorithms  for  enforcing 
nonnegativity  constraints  in  scientific  computation.  Special  emphasis  is  placed  on  such 
constraints  in  least  squares  computations  in  numerical  linear  algebra  and  in  nonlinear 
optimization.  Techniques  involving  nonnegative  low-rank  matrix  and  tensor  factor¬ 
izations  are  also  emphasized.  Details  are  provided  for  some  important  classical  and 
modern  applications  in  science  and  engineering.  For  completeness,  this  report  also  in¬ 
cludes  an  effort  toward  a  literature  survey  of  the  various  algorithms  and  applications  of 
nonnegativity  constraints  in  numerical  analysis.  Also  see  4  in  section  §2.1  and  reference 
13  in  section  §2.3. 

Tensor-Based  Analysis  for  Material  Identification  from  Hyperspectral  Data 

An  important  and  well  studied  problem  in  hyperspectral  image  data  applications  is  to 
identify  materials  present  in  the  object  or  scene  being  imaged  and  to  quantify  their 
abundance  in  the  mixture.  Due  to  the  increasing  quantity  of  data  usually  encountered 
in  hyperspectral  datasets,  effective  data  compression  is  also  an  important  consideration. 
In  references  2  of  section  §2.1  and  reference  11  of  section  §2.2,  we  develop  novel  methods 
based  on  tensor  analysis  that  focus  on  all  three  of  these  goals:  material  identification, 
material  abundance  estimation,  and  data  compression.  Test  results  are  reported  in  all 
three  perspectives. 
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4.2  Technical  Approach 
4.2.1  Pupil  Phase  Engineering 


Pupil  phase  engineering  (PPE),  as  it  appeared  in  Prasad,  et.  al.3  was  developed  as  a  con¬ 
strained  optimization  problem  to  minimize  sensitivity  of  the  PSF  to  nris-focus,  while  main¬ 
taining  a  restorability  metric  above  a  minimal  required  threshold.  A  Fisher  Information- 
based  metric  is  ideal  for  measuring  sensitivity  of  a  property  (e.g.,  the  PSF)  to  a  specific 
parameter,  (e.g.,  mis- focus).  Several  restorability  metrics  are  possible,  however  a  metric 
based  on  the  Strehl  ratio  proved  to  be  a  practical  choice. 

Pupil  phase  engineering  is  a  more  rational  approach  than  ad-hoc  choices  of  phase  mask,  since 
it  brings  numerical  constrained  optimization  techniques  to  bear  on  a  space  of  possible  phase 
masks,  arriving  at  an  optimum  solution  with  respect  to  a  well  defined  metric. 

For  purposes  of  extended  depth  of  focus  in  a  biometric  identification  system,  the  Fisher 
information  metric  is  minimized  within  a  constrained  set  of  phase  masks  ensuring  adequate 
restorability.  Minimizing  the  Fisher  Information  metric  minimizes  the  sensitivity  to  focusing 
errors. 

Let  u  =  (v,  w)  denote  an  image  plane  coordinate  vector,  and  let  r  =  (x,  y)  denote  a  pupil 
plane  coordinate  vector.  The  PSF  h( u)  in  the  image  plane  is  given  by: 

h(u)  =  |  AT(u)|2,  (1) 


where  K  is  the  following  pupil  integral  involving  the  total  pupil  phase,  [^(r)  +  0(r)]: 


K{  u) 


1 

A  fVA 


J  drP(r)ei[^u'r+^r)+e{r)] . 


(2) 


In  this  expression,  P( r)  is  the  pupil  function,  equal  to  1  inside  the  pupil  and  0  outside  for  a 
clear  pupil,  u-r  denotes  the  dot  product  of  vectors  u  and  r,  A  is  the  wavelength  of  illuminating 
light,  /  is  the  focal  length,  and  A  is  the  area  of  the  pupil.  For  the  purpose  of  extended  depth 
of  focus  for  biometric  identification  systems,  the  uncompensated  pupil  phase  0(r)  is  of  form 
t(x2  +  y 2). 

Pupil  phase  engineering  seeks  a  phase  mask  that  allows  the  greatest  possible  insensitivity 
of  the  phase-encoded  optical  image  to  focus  variation  without  unacceptably  compromising 
the  digital  restorability  of  that  image.  A  Fisher  information-based  metric  measuring  the 
sensitivity  of  the  PSF  to  defocus  may  be  defined  as: 


Ji(r) 


1 


h(  u,  r) 


<91n/r(u,  t)12 


du 


(3) 


Jj(t)  is  a  weighted  average  of  the  square  (in  log  scale)  of  the  rate  of  change  of  the  point 
spread  function  with  respect  to  the  defocus  parameter  r.  If  Jj{t)  is  small,  we  expect  the 

3S.  Prasad,  V.P.  Pauca,  R.J.  Plemmons,  T.C.  Torgersen,  and  J.  van  der  Gracht.  Pupil-phase  optimiza¬ 
tion  for  extended-focus  aberration-corrected  imaging  systems.  In  Proc.  SPIE,  Advanced  Signal  Processing 
Algorithms,  Architectures,  and  Implementations  XIV,  volume  5559,  pages  335-345,  Denver,  CO,  Aug.  2004. 
SPIE,  SPIE. 
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PSF  to  change  very  little  in  a  neighborhood  of  r.  A  natural  choice  of  optimization  metric  is 
given  by  the  following  integral  of  [ J/(r )]2: 

I(r0)=  r  [Ji(t)}2  dr,  (4) 

J  —TO 

where  (— to,  to)  denotes  the  (symmetric)  range  of  defocus  parameters  of  interest.  The  choice 
of  To  represents  an  important  opportunity  to  integrate  the  mission  requirements  (required 
depth  of  focus)  of  the  imaging  task  directly  into  the  optimization  step. 

Minimizing  I(tq)  also  admits  phase  masks  for  which  the  corresponding  restoration  problem  is 
too  ill-conditioned  for  practical  use.  A  penalty,  or  barrier  function  is  utilized  to  constrain  the 
minimization  of  / (tq)  to  allow  only  phase  masks  whose  corresponding  restoration  problems 
are  sufficiently  well  conditioned. 


4.2.2  Phase  Encoding  to  Mitigate  a  Pulsed  Laser  Attack  on  a  CCD  Camera 

The  problem  assumes  a  frequency- agile,  pulsed  laser  capable  adversary  with  intent  on  dam¬ 
aging  or  disabling  a  CCD  surveillance  camera.  One  possible  counter-measure  is  to  use  phase 
encoding  to  spread  the  incoming  beam  to  reduce  the  destructive  energy  delivered  to  any 
one  point  on  the  detector  surface.  Image  reconstruction  can  (in  theory)  be  then  used  to 
reconstruct  an  image  suitable  for  the  surveillance  application  of  interest.  A  number  of  phase 
element  design  strategies  have  been  considered  including  traditional  cubic  phase-mask,  piece- 
wise  linear,  pseudo-random  phase  ,  and  modified  pupil  phase  engineering.  It  is  straightfor¬ 
ward  to  design  a  piece-wise  linear  phase  element  which  takes  an  incoming  point  source  and 
spreads  it  into  100  equally  spaced  points  containing  equal  energy;  this  guarantees  the  required 
Strehl  ratio  of  0.01.  Figure  1  (left)  illustrates  a  point  spread  function  using  this  approach. 
A  trade-off  exists  between  held  of  view  and  quality  of  restoration.  If  the  held  of  view  is 
appropriately  masked,  sub-images  spread  across  sub-regions  of  the  detector  do  not  overlap, 
and  may  be  treated  as  independently  sampled  images  of  the  same  scene.  The  generalized 
sampling  theorem  (see  reference  11)  can  then  be  used  to  re-construct  a  high-quality  image. 
A  simulation  study  has  shown  this  approach  offers  very  good  restoration  in  the  presence  of 
noise,  but  sacrihces  held  of  view  to  attain  the  desired  mitigation. 
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Figure  1:  Point  spread  function  (left)  and  corresponding  MTF  (right)  for  piece-wise  linear 
phase  perturbation 

A  summary  of  our  results  for  this  problem  can  be  found  in  results  [6]. 
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4.2.3  Lenslet  Array  Imaging 


Image  formation  using  an  array  of  lenslets  focused  on  a  single  detector  provides  the  raw 
data  to  reconstruct  a  high  resolution  image  using  the  generalized  sampling  theorem.  In  this 
context,  it  is  important  to  clarify  the  distinction  between  digital  super-resolution  and  optical 
super-resolution.  Digital  super-resolution  extends  the  image  resolution  beyond  the  sampling 
rate  of  the  detector.  Given  multiple  frames  of  the  same  scene,  recorded  with  appropriate 
sub-pixel  shifts  among  the  frames,  a  well-understood  theory  shows  that  a  high  resolution 
image  can  be  computed.  Optical  super-resolution  extends  the  image  resolution  beyond  the 
diffraction  limit  of  the  optical  system.  Recent  information  theoretic  results  show  that  optical 
super-resolution  is  not  possible  without  additional  sources  of  information,  such  as  support 
constraints,  non-negativity  constraints,  etc. 

Lenslet  Camera  Image  Formation  Model 

til  •  • 

A  low-resolution  image,  gj  formed  by  the  j  lenslet  is  given  by: 

9j  =  DHjSjf  +  Vj,  (5) 


where, 


•  /  is  a  target  high-resolution  representation  for  the  object, 

•  Sj  represents  the  translation  of  image  /  due  to  the  relative  position  of  lenslet  j  with 
respect  to  the  object  and  the  detector, 

•  Hj  is  a  blurring  operator  associated  with  lenslet  j, 

•  D  is  a  decimation  operator  which  represents  the  transformation  from  the  target  high 
resolution  to  the  (low)  resolution  of  the  detector, 

•  gj  is  the  low-resolution  image  associated  with  lenslet  j  and, 

iL 

•  rjj  describes  a  noise  process  associated  with  the  j  J  image. 

Inverse  Problem  Definition 

For  an  n-lenslet  camera  we  form  a  least-squares  functional: 


J(f) 


DRiSi  ' 

9i 

dh2s2 

92 

_  DHnSn  _ 

9n 

whose  minimum  gives  the  desired  estimate  /  for  the  high 
resolution  image  /.  I.e., 


(6) 


/  =  argminj  {«/(/)} 


(7) 


Inverse  Problem  Structure 

When  a  pixel  in  /  is  shifted,  it  can  only  overlap  a  small  number  of  neighbors.  This  implies 
there  are  only  a  small  number  of  non-zero  elements  in  each  column  of  Sj.  Conservation  of 
light  flux  implies  Sj  must  be  column  stochastic  except  at  in  columns  corresponding  to  the 
boundary  of  the  detector.  Hj  is  circulant,  and  therefore  efficient  Fourier  transforms  can  be 
used  to  evaluate  the  forward  model.  The  down-sample  operator  D  only  integrates  over  a 
small  number  neighboring  pixels  in  HjSjf .  This  implies  there  are  only  a  small  number  of 
non-zero,  elements  in  each  row  of  D.  The  structure  of  the  inverse  problem  implies  that  the 
forward  model  should  be  kept  in  factored  form  and  that  an  iterative  approach  is  most 
favorable.  The  current  implementation  described  in  [12]  uses  the  well-known  CGLS  method. 

Application  to  Biometric  Identification 

Figure  (2)  compares  the  performance  of  a  single  lenslet  with  our  digital  restoration  using 
multiple  lenslets  under  two  distinct  processing  options.  The  horizontal  axis  represents  the 
number  of  iterations  used  in  the  restoration  method.  The  vertical  axis  represents  the 
Hamming  distance  between  the  Daugman  encodings  of  two  iris  image  captures:  a 
high-resolution  reference  image,  and  an  image  syntheized  from  data  produced  by  a 
prototype  camera.  The  dashed  horizontal  line  at  0.32  on  the  vertical  axis  represents  the 
threshold  for  successful  biometric  identification  based  on  an  iris  match.  The  graph  clearly 
show  a  favorable  result  when  selecting  a  subset  comprised  of  the  best  sub-images. 
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Figure  2:  Comparison  of  single  lenslet  performance  with  multiple  lenslet  performance. 


9 


4.2.4  Non-negative  Matrix  Factorizaton  and  Non-negative  Tensor 
Factorization 

Three  major  objectives  in  processing  hyperspectral  image  data  of  an  object  (target)  are: 
data  compression,  spectral  signature  identification  of  constituent  materials,  and 
determination  of  their  corresponding  fractional  abundances.  In  reference  [11],  the  authors 
propose  a  novel  approach  to  processing  hyperspectral  data  using  non-negative  tensor 
factorization  (NTF),  which  reduces  a  large  tensor  into  three  factor  matrices;  the  Khatri-Rao 
product  of  these  factors  approximates  the  original  tensor.  This  approach  preserves  physical 
characteristics  of  the  data  such  as  nonnegativity  and  can  be  used  to  satisfy  all  three  major 
objectives.  Test  results  are  reported  in  reference  [11]  for  space  object  identification. 

Nonnegative  Factorizations 

In  Nonnegative  Matrix  Factorization  (NMF)  and  m  x  n  (nonnegative)  mixed  data  matrix 
X  is  approximately  factored  into  a  product  of  two  nonnegative  rank-A;  matrices,  with  k 
small  compared  to  m  and  n,  X  ~WH .  This  factorization  has  the  advantage  that  W  and  H 
can  provide  a  physically  realizable  representation  of  the  mixed  data.  NMF  is  widely  used  in 
a  variety  of  applications,  including  air  quality  control,  image  and  spectral  data  processing, 
text  mining,  chemometric  analysis,  neural  learning  processes,  sound  recognition,  remote 
sensing,  and  object  characterization.  Nonnegative  Tensor  Factorization  (NTF)  is  a  natural 
extension  of  NMF  to  higher  dimensional  data.  In  NTF,  high-dimensional  data,  such  as 
hyperspectral  or  other  image  cubes  is  approximated  by  a  sum  of  rank  one  nonnegative 
tensors.  The  algorithms  given  below  combine  features  of  both  NMF  and  NTF  methods. 


Definition  1:  We  define  a  non-negative  rank-/c  decomposition  of  the  tensor  r  as: 


nun 

x(t)  ,y(i)  ,2(0 


r  — 


Y  X (i)  o  y W  o 


i—  1 


subject  to 

x(i)  >  0  ,  yW  >  o  ,  Z(i)  >  0 

where  r  €  nD'xD2xDs, x(i)  €  ■RD1  ?  y(i)  G  nD2  ?  and  z(i)  €  nD3 


A  graphical  depiction  of  the  decomposition  problem  is  given  in  figure  3  below. 
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Figure  3: 
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NTF  Algorithm 


•  Group  Xj’s,  yi  s,  and  Zi  s  as  columns  in  X  €  lZ^lXk,  Y  €  TZ+2Xk,  and  Z  €  lZ+3Xk 
respectively. 

•  Initialize  A  and  V  by  Nonnegative  Matrix  Factorization  of  the  mean  slice, 

min  II  A  —  XYT  \  \  \ 

XY  11  llf 

where  A  is  the  mean  of  r  across  the  third  dimension. 

•  Iterative  Tri- Alternating  Minimization 

1.  Fix  t.  X .  Y,  and  fit  Z  by  solving  a  NMF  problem  using  a  projected  gradient 
descent  algorithm. 

min  1 1  Tz  —  CzZT  ||p 
z  11 

where  Cz  =  XQ  Y  g  KDlD2Xk,  Tz  is  the  unfolding  tensor  across  the  third 
dimension,  and  Q  represents  the  Khatri-Rao  product. 

2.  Fix  r,  X,  Z,  and  fit  Y  as  above. 

3.  Fix  r,  Y,  Z ,  and  fit  Z  as  above. 

The  algorithm  given  above  has  been  shown  in  simulation  to  effectively  identify  the 
fractional  abundances  present  in  a  computer-generated  data  model  of  the  Hubble  Space 
Telescope.  A  data  volume  of  177  x  193  x  100  was  successfully  processed. 


5  Technology  Transfer 

5.1  Mitigating  Pulsed  Laser  Attack  on  CCD-based  Cameras 

The  WFU  imaging  group  was  contacted  in  April  2008  by  Agiltron  Corporation  to 
cooperatively  develop  an  SBIR  proposal  entitled  “Utilizing  Computational  Imaging  for 
Laser  Intensity  Reduction  at  CCD  Focal  Planes”4.  This  collaboration  facilitates  technology 
transfer  and  commercialization  of  Pupil  Phase  Engineering  (PPE)  methods  developed  by 
the  WFU  imaging  group  and  Dr.  Sudhakar  Prasad  of  the  University  of  New  Mexico.  The 
collaborators  wish  to  thank  the  Army  Research  Office  (ARO)  for  their  generous  support 
during  the  development  of  PPE  methods. 

Todd  Torgersen  and  Joe  van  der  Gracht  of  Holospex,  Inc.  assisted  Agiltron  with  the 
technical  portion  of  the  phase  I  proposal,  and  served  as  consultants  during  the  performance 
period.  Agiltron  and  WFU  were  invited  to  submit  a  phase  II  proposal,  and  submitted  one 
in  September  2009.  While  this  proposal  was  not  selected  for  phase  two  funding,  the 
investigators  showed  that  their  approach  is  viable. 

The  problem  of  interest  is  to  apply  computational  imaging5  to  mitigate  the  effect  of  a  pulse 
laser  attack  on  a  CCD  based  camera.  The  basic  approach  is  to  choose  a  pupil-phase 

4See  topic  A08-064,  URL:  http://www.armysbir.com/awards/sbir_08ph.l_company.htm 
5  Optical  design  combined  with  digital  post-processing. 
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perturbation  which  will  spread  out  the  incoming  energy  such  that  the  PSF  peak  energy 
impacting  any  single  detector  pixel  is  reduced  by  a  factor  of  100.  The  assumption  of  a 
frequency-agile,  pulsed  laser  capable  adversary  excludes  solutions  based  on  simple 
band-pass  filtering,  electro-mechanical  shutters,  or  even  digitally  controlled  transmissive 
LCD  elements. 

Several  pupil  phase  encoding  elements  were  investigated,  including  a  traditional  cubic 
phase  pertubation,  a  piece-wise  linear  element  relating  to  sub-aperture  techniques,  and  a 
pseudo-random  based  phase  pertubation.  These  approaches  are  discussed  in  [6]. 

5.2  Rapid  Assessment  of  Thermal  Injury 

The  co-investigators  are  currently  collaborating  with  Dr.  J.  Molner  at  the  Wake  Forest 
University  School  of  Medicine  to  adapt  multi-spectral  lenslet  array  imaging  techniques  to 
the  problem  of  rapid  assessment  of  thermal  injuries.  Unfortunately,  thermal  injuries  are 
often  progressively  widening  (or  deepening).  Injured  (but  still  living)  tissue  near  “ground 
zero”  may  receive  reduced  blood  flow  in  the  days  following  the  injury,  leading  to  an 
expanding  area  of  non-viable  tissue.  Traditional  treatments  may  require  hospitalization, 
heavy  medication,  and  very  low  patient  mobility  for  several  days  before  the  correct  surgical 
(e.g.,  skin  graft)  boundaries  become  apparent.  The  goal  of  this  study  is  to  evaluate  several 
possible  technologies  which  may  predict  viable  and  non-viable  regions  at  the  earliest 
possible  time.  Multi-spectral  imaging  allows  us  (in  theory)  to  map  regions  of  reduced  (or 
increased)  blood  flow  based  on  localized  oxygenation  levels.  The  spectral  response  of 
hemoglobin  and  oxy-hemoglobin  differ  markedly  at  several  wavelengths,  thus  the  spectral 
response  is  well  correlated  to  percentage  of  blood  oxygenation.  The  challenges  posed  by 
this  work  include  1)  highly  accurate  calibration  of  all  aspects  of  data  acquisition,  2)  the 
effect  of  skin  pigmentation,  and  3)  integration  of  accurate  scattering  models.  The  results  of 
this  investigation  are  found  in  “Assessment  of  hemodynamic  changes  of  pig  skin  in  the 
post-burn  period  using  a  multi-spectral  multi-aperture  camera”,  WFU  technical  report. 
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